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Abstract 

Precision Satellite to Satellite Tracking can be used 
to obtain high precision and resolution maps of the geoid. 
A method is demonstrated to use data in a limited region to 
map the geopotential at the satellite altitude. An inverse 
method is used to downward continue the potential to the 
earths surface. The method is designed for both satellites 
in the same low orbit. 

Introduction 

With the general acceptance of the need for improved 
global geoid and gravity data for both Solid Earth 
Geophysics and Physical Oceanography investigations, the 
fact that this data can only be provided with a close earth 
satellite measurement system, and the realization that at 
least one satellite measurement system is now able to 
provide observations with the required accuracy (NAS/NRC 
1979 report) , attention must be focused on the details of 
extracting the geoid and gravity data from these 
observations in an efficient way that will also be suitable 
for these investigations. The view taken here is that the 
Gravity Research Mission (GRM) , a Satellite to Satellite 
Tracking Mission (SST) will satisfy these objectives. 
Generally we assume that the SST mission will have two 
surface force compensated satellites in similar low altitude 
orbits (160 to 180 km) separated by 300 to 1000 km and a 
tracking system with an accuracy of 1 micron/sec (l///sec = 
10**-6 m/sec) for 10 second averages. This analysis 
concentrates on this generic mission. 

The concept of the SST mission is simply to calculate 
the gravitatinal force acting on a spacecraft from changes 
in its measured velocity. The satellite itself is the 
sensor, and its velocity, the observable. This is pictured 
simplistically as follows. Consider a satellite in orbit 
around the earth approaching a region of excess mass. As 
the satellite approaches, it is accelerated toward the mass, 
and after passing the mass, it is decelerated. By measuring 
the time history of the velocity variation, an estimate of 
the magnitude and position of the mass excess can be 
deduced. Of course the actual situation is much more 
complex, for several reasons: the structure of the earth's 
mass distribution is very complicated, other forces act on 
the satellite, only one component of the satellite velocity 
is measured from another satellite, and the observations 
contain errors. In the SST concept, the second satellite 
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could be very high, say in geosynchronous orbit frhe 
high-low configuration). Alternatively, the second 
satellite could be in nearly the same orbit trailing the 
first, low satellite (the low-low configuration) . In the 
low-low case, the second satellite would experience similar 
velocity changes, but at a later time. 

SST tracking has been realized in a number of missions. 
The earliest example of doppler tracking a satellite from a 
point not on the body being studied is the mapping of the 
lunar gravity field by means of lunar orbiters tracked from 
the earth (Muller and Sjogren, 1968) . Since this remarkable 
success, there have been notable doppler tracking 
experiments for earth satellites: the tracking of Geos-3 
from ATS- 6 geosynchronous satellite, the tracking of the 
Apollo Command Service Module (CSM) from the ATS-6, and the 
tracking of the Apollo CSM from a small satellite trailing 
it in the same orbit. The first two are examples of the 
high-low configuration, and the last is an example of the 
low-low configuration. 

The classical method of determination of geopotential 
information using perturbation methods, and spherical 
harmonic representation could be employed (Gaposchkin, 
Kaula, Guire & Newton, Balmino, Lerch et al etc). However, 
the increased resolution makes this approach unattractive. 
For example, assuming a desired resolution of 1 degree, i.e. 
a spherical harmonic expansion to degree and order 360, 
requires the simultaneous recovery of almost 130000 
coefficients, in addition to the orbital parameters. Direct 
recovery of gravity anomalies from SST data has also been 
investigated (eg Hajela 1977) . Some view the problem as an 
estimation problem in a finite-dimension vector space of 
parameters (Vonbun et al., 1978; Schwarz, 1970). Others 
have treated the problem as mapping the acceleration field 
at the satellites altitude (Muller and Sjorgren, 1968; Marsh 
etal. 1977). Studies have also been done on mapping the 
acceleration or potential at the earth's surface (Rummel, 
1975, Rummel et al, 1978). Some consideration has also been 
given to use of fourier techniques (Colombo, Kaula) , and 
they show promise of significant economies in computer 
resources. The development presented here takes the view, 
that we can use the data directly to recover the desired 
quantity. It is analgous to an approach used by Muller & 
Sjogren 1968, and Sjogren 1976. 

We propose that the SST mission be viewed as a mapping 
mission. That is, the result will be maps of the geoid or 
gravity, as contrasted to determination of spherical 
harmonic or fourier coefficients. Of course (digital) maps 
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can be used to derive these coefficients. As a mapping 
mission, the computation can be done area by area. This 
avoids the enormous computer costs of simultaneously 
recovering all the geopotential information and the 
trajectories. 

The process of mapping the SST signal to a geoid would 
have a number of steps: 

1) With a chosen geopotential model, compute the range 
rate residuals from the best trajectory available. This 
need only be done once. 

2) With a suitable algorithm, map the differential 
velocity ( £v) to the anomalous potential ( AT) at the 
satellite altitude. 

3) Regularize the potential along the satellite tracks 
to the potential on a sphere near the satellite altitude. 

4) Downward continue the potential from the reference 
sphere at the satellite altitude to the ellipsoid. 

5) Correct the SST signal with the (revised) anomalous 
potential, and repeat from step 2. 

The approach has a number of advantages. First, it 
separates the inferance of the geopotential from the 
downward continuation problem. One can evaluate the 
geopotential at satellite altitude independently from its 
use on the ellipsoid. For some uses, e.g. precision orbit 
computation, the result can be used directly. Second, the 
mapping can be done locally. One can develop a map of the 
geopotential for a region, without solving the complete 
problem. This results in more manageable computations, and 
control on the ingrediants of the computation. Third it 
does not require costly recalculation of trajectories using 
the improved potential at each iteration. It allows 
proceeding arc by arc, and comparison of the restults at 
each step with other data. 

What follows is the formulation of an analytical 
solution of the problem of calculating satellite 
perturbations caused by the anomalous potential. Such a 
solution provides insight into the proposed measurement with 
significant computational economies. It gives a direct way 
to calculate analytically partial derivatives of the 
observable (velocity) the respect to the desired end product 
(the geoid) . The analytical solution also allows the kernel 
function in the generalized inversion to be calculated 
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directly. Furthermore we have investigated the use* of 
generalized inverse (Backus and Gilbert, 1967,1968, 1970; 
Parker, 1975) to convert the observed velocity measurements 
to a mapping of surface values of geoid height. The 
ill-posed boundary value problem and the unstable downward 
continuation problem are addressed in a way that the 
additional assumptions used to obtain the solutions are 
clearly identified and an error estimate can be found. 
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The Range Rate Equation 

We can begin with a geometrical discussion of the range 
equation. In general the distance between two satellites 
( A ) with positions r t * , <f>. , , i=l f 2 can be written 


2 2 2 

A = r x + r 2 -2r^ r 2 cos if/ 


(1 


where the (geocentric) angle \j/ between them is given by 


cos 0 = sin^ x sin^ + cos <^cos ^cos ( X 2 ) 


(2 


The range rate equation is easily found to be 


0 (y\ - -V if/ 2 *% » jyysWt ^ 


(3 


Equations 1,2, and 3 are valid for any satellite positions, 
and are not restricted to the low-low configuration. We 
proceed by assuming a reference orbit (r 0 ; ,<&>;, X*:, i=l,2) and 
potential U and determine corrections (ir t * r6<f>i 1 6 Xi , i=l , 2) 
to the orbit. The dynamical equations enter when relating 
these corrections to corrections to the potential ( T) . 

When dealing with perturbations of a reference orbit we have 

A = A 0 ♦- 6A 

A = *3A < 4 

+1 = 

Considering 6 A ,6<\ to be quantities of first 

order we have 


A = 


( f » « ~ K ) rt t( ^ (5 
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- ) St»* 

+( Sr*4 f *”*4x Or* 0*»“3O} ^ +%, 


This range rate equation is valid for any combination of 
high and low satellites. The perturbations 6r,6<t> , and 6X 
are just those calculated by the theory given in the next 
section when it is applies in a reference frame defined by 
the earth's equator and prime meridian. Therefore this 
theory is applicable to any combination of satellites. If 
we do make the assumptions of two low small eccentricity 
satellites, then we obtain the usual relation (r* = R + dr; 
,i=l,2) 


t * <5t? 

o A = — — — — 

A 


(7 


Proceeding to treat the two coplanar low satellite 
case, it is convenient to use the perturbation solution in a 
reference frame defined by the mean orbital plane of the two 
satellites and its equator crossing. In this case the 
perturbations in radius are the same, and the separation 
angle becomes the difference in true anomaly of the two 
satellites which is the same as the difference in longitude 
[66) of the theory. In this case 
( 4>- =0, i=l ,2) , ip = X- A = 6> - 6^ «1 and we have 


6 A = 6Xj = R(d<j- 66 x ) (8 


The solution (eq 33) specialized for small eccentricity with 
R=a is 


T 


Zr> 6r( T) 



66 = 


*na 


(9 


K 
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where dr(T) is given in equation 27. We now have the range 
rate equation 

T{.? TM i \ 

6 A = - " 2n( 6'CLO - dtix) J (10 


which becomes our observation equation. It will be shown 
that the term n jc is smaller than T/na. Then the /irst 
terms provide the desired mapping of the observable ( 6 A ) to 
the potential difference T(l)-T(2). Now dr is a rather 
complicated integral function of T, depending on the past 
motion of the satellite. It is critical to recognize that 
it is the complex dependence of Sr on T that mandantes an 
iterative proceedure to obtain T. This arises due to the 
fact that dT/ dr occurs in the expression for Sr. 

To summarize the discussion to this point, we have 
presented an iterative scheme that uses residuals in 
satellite to satellite range rate ( sA ) with respect to a 
reference trajectory and geopotential model (Uo) . The 
scheme proceeds in two stages. A number of crossing tracks 
are regularized by imposing that the crossing points, 
analytically continued to the same reference surface (a 
sphere) , have the same value of potential. This computation 
is analagous to crossing arc analysis used with satellite to 
sea surface altimeter data. This regularization cannot be 
done perfectly due to the 5r terms in equation 9, which 
require a_priori knowledge of the potential to be computed. 
Therefore the iteration proceeds by using the potential 
(T„-| ) calculated in iteration n-1 to calculate 5r and 
d%/ dr for iteration n. The resulting potential (T) can 
then be analytically continued to another surface using 
methods based on inverse theory. 
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Theory of Satellite Motion 

In discussing the recovery of gravity anomalies ( g) 
or geoid heights (N) with the SST approach, it is 
convenient, if possible, to have analytical formulae 
relating satellite position and velocity (the state vector) 
to the desired quantity. Such formulae also enable the 
sensitivity, or partial derivatives, of gravity anomalies 
and geoid heights to be more easily obtained; otherwise one 
would have to resort to costly numerical methods, with 
consequent loss of generality and insight. 

For discussion purposes, we start with the theorem of 
conservation of energy of the satellite orbit. Of course, 
drag and radiation pressure perturbations would also have to 
be taken into account if the following relationships were 
used for actual data. One can obtain the same result for 
circular orbits by using the conservation of angular 
momentum. Let the total potential be represented by , 
which, for convenience, can be separated into a reference 
potential u and an anomalous potential T, i.e. 


^ = U + T 


(11 


In standard physical geodesy notation, U is the normal 
potential corresponding to that of a reference ellipsoid; 
however here we prefer to view it as a reference or adopted 
potential, with T being the remaining (unmodeled) part that 
we seek. Some confusion is bound to occur because the 
anomalous potential, or disturbing function, is generally 
denoted by R in celestial mechanics and by T in Physical 
Geodesy. If we write the kinetic energy as mv i /2, with 
the vector component of velocity along track (v u ) , cross 
track (v ) and radial (v r ) , then 


t ( v w + + v r ) + Sf' = constant 


(12 


For a satellite with small eccentricity, we can treat the 
along track velocity and the unperturbed velocity vo, which 
gives v u = vo + $ v u r %, = Sv w , and v r = 8v r . 
Therefore, to first order in small quantities. 
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5v u = T/vo (13 


This formula was first derived by Wolf (1969) and 
subsequently used by many others. For small eccentricity 
satellites, equation (13) gives the change in along track 
velocity as a function of the potential in space (i.e. 
position and time) to about 10%, the errors arising from the 
change in radial distance owing to T, interacting with 0. 

Equation 13, a direct mapping of the potential into the 
along track velocity, could be used as a first approximation 
for inverting observations of 5v w to determine T. 
However, significant perturbations in 8v w and 6v r also 
contain information about the disturbing potential T, so to 
obtain general formulae for the three components, we must 
solve the dynamical equations. This is done by using 
established proceedures in planetary theory for special 
perturbations. 

Following the treatment in Brouwer and Clemence (1961) , 
we can write the equations of motion as 


clt X 


X 

■ s 

£T 

*> 


cl\ 

+ >*- 


OT 

(14 

cL*^ 


± _ 
T* 

52T 

St 


v - 

%o 3 

GM = 3.986x10 cm 

sec ^ and 


z 

r 

II 

% 

+ 

* A 2. 
y + z 


(15 


For the moment, we assume T excludes only the central force 
term from . If the coordinates are r, 6 , and <t> , 
where 0 is the longitude from the reference (x y plane) 
equator crossing in the reference plane, and <f> is the 
latitude, then the equivalent differential equations of 
motion are 
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(16 


Now these equations are most commonly used with respect to 
the earths equator. However, they are valid for any 
inertial frame. For our purposes it is convenient to adopt 
the orbital plane, at the beginning of the track, as the 
reference equator. In this frame, 0 becomes the along 
track coordinate and 0 the cross track coordinate. The 
origin, zero point, of course is arbitrary, and we chose the 
beginning of the data for convenience. We now must express 
the potential in this orbital frame. 


To obtain an expression for the potential (T) , and its 
gradiants, as a function of position (and time), we begin by 
assuming the potential on the earth's surface at r=R. Since 
T(R, <{> , X ) given on a sphere is the Dirichlet boundary 
condition for Laplace's equation external to the sphere r=R, 
we can use the basic results from potential theory to obtain 
T at any point in space. Therefore, we can expand T in 
terms of orthonormal base functions (associated Legendre 
functions) as 


T(R ,0rX) 



C cos mX + 



sin mX ) Pj^sin <f> ) 


(17 


which, using the properties of solutions of Laplace's 
equation in spherical coordinaes, can be upward continued as 


T(r,0,X) 



C^cos mx + 


S. sin mX 

Ih, 


)P. (sin 4 > ) 


(18 




Alternatively, since T is harmonic in space, we can use 
Poisson's integral to obtain 



1 f 
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T(r ,<f>,k) 



TIAz£lD ^ JL<t>‘ JL>! 


(19 


Wa. 


where p is 
Rf <t>' f X and 
written as 


the distance between the integration point 
the sample point r , <f> , \ ; it is often 


- r*“ + R*” - 2Rr cos ip 


(20 


ip being the central angle, which can be expressed as 


cosip = sin0 sin0 / + cos<£ cos0 / cos ( X - A 7 ) 


(21 


We recall that the relation between the anomalous 
potential T and the geoid height N by Brun's formula is 


N = T/y 


(22 


where y is the normal gravity on the reference surface. 
Furthermore, the relation between gravity anomalies Aq and 
T is given by 


Aq 


£T _Z T 
ay* y* 


(23 


and we can formally obtain N via Stokes' equation 


4-7T y J 


S ( ip ) Aq da 


N 


(24 
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and Molodenskii ' s equation 



_L_ f 

lfc 7 T R 5 J 


JLilLi jlo 


(25 


Both equations (18) and (19) are suitable for use with 
the equations of motion (16) . The spherical harmonic form 

(18) is equivalent to classical perturbation analysis and 
thus is not of particular interest here. However, equation 

(19) will allow mapping of the velocity directly to a 
surface potential. 


Following the approach in Brouwer and Clemence (1961) 
we can write the perturbations in r ,6 , and <f> , as 
follows. Let 


X = 


-C 


H 


cLT 


£_r 

d <* 


] 


Y 



(26 


1 


•C 


where p=a(l - e 1 * ) , a is the satellites semi major axis, e 
the eccentricity, and r corresponds to the unperturbed 
motion. We have 


6 * * y X 
60 = fl 


In equation (26) we interpret 
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K ' 


dT = 


dT 

ar 

dT 

S’*- 


dL-r +■ 
•l-x. +- 


d T 

do' 

dT 

0 


<LQ + 
«L^ + 


dT 

d<f> 

dT 

d* 


•«40 


= (VT . Tf )o4t 


(28 


Therefore, for example the integral in the first part of 
equation (27) could be obtained by expressing r and T in 
terms of the true anomaly f, and it becomes 



(29 


Other choices for the variable of integration are possible; 
for example, we could use the mean anomaly M, or the 
eccentric anomaly E. In equation (27) we have written them 
as integrals in f, although the same choices obtain. The 
form, using the true anomaly (f) is the one implemented in 
the following analysis. Note that the notation "(f)", due 
to Hansen, means that the variable is held fixed in the 
integration but takes its proper value when the perturbation 
is evaluated. These expressions then have the mathematical 
form of a convolution. 


Equations (26) and (27) are quite general and depend on 
finding a once differentiable potential function T. The 
observable velocity can, of course, be found from the 
appropriate linear combination of derivatives of equation 
(27) : 


-7T 

r 


-xr 

u 


- i? L /x-tw-fl" 

"i? • f '/ T " • 


V 

W cLt 


(30 


n. 







As a mater of interest, these expressions can be simplified 
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for very small eccentricity: 

X ' i [^T *-#] 


Y - \ d JL 

1 “ ^ V 2 * 5 M 


2 = 


1 d T 


'h’oJ d St- 


and 


r 


■r -n. 


/ x «•»£ (m) - H] 4.1^1 


(31 


< *U* rr *7L 

U 


■A 


dM — 2. 5v* 


"'ll/ - A <«n£ CM') — 


(32 


We can rederive equation (13) from these equations. 


"XT r t\ 


/ T - 7. % -r) s T* — Xxl 

\ >V ' ^ 


(33 


where we can make the identification vo = na. The second 
term can be written: 


~yi%T = - £22 s* ^ 

.3 


= -4- f— )l«v 

"V; a.’- /I 


rc«x (34 
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which is simply the change in velocity due to the 
perturbation 5r interacting with the central force term, as 
was pointed out above. 

This development can be employed in a number of ways. 
The particular computer implementation made the following 
choices. The perturbations were numerically computed using 
the Poisson integral (19) form of the potential. The 
reference trajectory was given and the gradiants were 
evaluated for each data point. Then a simple Simpsons Rule 
quadrature was used to calculate the integrals for each data 
point. Of course, the inital values of the integral were 
assumed zero. Any constant of integration is applied during 
the regularization, when all the track crossings are 
adjusted, and to have zero mean. The gradiants were 
computed in an earth fixed system, and then transformed to 
the adopted orbital plane, along track, cross track, and 
radial as follows. In the terrestrial frame we have: 


£1= JL /[ - -t- * (<V» A T JLa' (35 

ar An J c * J p* 


The gradiant vector can be expressed in the inertial frame 
with 

faT/r vn(j)d\ 1 


F 




a \ /'Cd<}> 


(36 


aT /at* 


where Rj(x), R a (x) , and R 3 (x) are three dimensional 
rotations about the x,y and z axes respectively, and $ is 
the siderial angle for the time in question (Gaposchkin 
1973) . The along track, cross track and radial components 
are given by 



If 
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(37 


where the along track, radial and cross track unit vectors 
are 


*« = v / W 

% = T/1r\ (38 

U u> = u r x V 1 u r X u ul 


The last important detail is the implementation of the 
Poisson integral. It is formally defined as an integral 
over the sphere. We assume that only the nearby area make 
nonzero contributions to the integral. This is based on the 
notion that since we are using a relatively high degree and 
order reference potential, (e.g. l=m=12) distant zones 
formally average to zero. Therefore, to limit the 
quadrature, the Poisson integral only includes surface 
elements less than 15 times the minimum distance from the 
evaluation point and the reference sphere, i.e. p < 
15(r-R) . 
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Downward Continuation and Inverse Theory. 

The two issues concerning the fundamental process of 
converting potential measurements at satellite altitudes to 
potential, i.e. geoid heights, at the earth's surface are 
the solution of a boundary value problem in potential theory 
with incomplete data on an undefined surface, and the 
downward continuation of the potential to the earth's 
surface in a manner that keeps the errors within bounds. 
The approach analysed here address both problems in one 
step. We will use generalized inverse methods (Bakus and 
Gilbert, 1967, 1968, 1970; Parker, 1977). This approach 
provides a means to incorporate error analysis, and to 
provide an arbitrary amount of smoothing. 

Briefly, observations of 6v are averages over some 
integration time, a property of the' instrumentation. This 
is typically on the order of 10 sec. These observations are 
given on a smooth curve, but the curve is not a suitable 
boundary for solution of Laplace's equation. Therefore, we 
must somehow express the observed 8v as if if were on a 
boundary. Also, the boundary values are given at a finite 
number of points, immediately precluding a complete 
solution. The limited number of points requires some 
additional assumptions about the potential field, either in 
terms of smoothness or as a_priori information. In 
addition, the observations will leave errors. The 
analytical continuation of the potential field is a 
stability issue (Bullard and Cooper, 1948): Any short 
wavelength error in the potential at satellite altitudes 
will grow arbitrarily large when downward continued, and the 
chosen method must control this error build up in a well 
defined way. 

One common practice is to limit the vector space of the 
representation by performing a least squares solution 
(Vonbun et al. 1978) . This has the advantage of a definite 
computation but does not identify the errors in the 
recovered anomalies as a function of wavenumber. 


In the study performed, the following outline of 
spectral expansion inverse theory, we follow Parker (1977) . 
Consider that we can write T at a point P* as 


/ 




T( » 


dor 


A 


T(P') 


(39 
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This is exactly the form of the Poisson integral. In 
equation (39) Pi is the point where the observation is made, 
and P^ are on the desired surface, which can be taken as the 
earth's surface, and the element of integration d is on 
that surface; G(P ,P )is now a kernel function. We can 
treat the observation errors by using 


Tj| = T/<t : 


(40 


where oi is the standard error of the observation T. The 
matrix with elements is defined as 

f 

I G(P. ,P^)G(P. , P^) da (41 

For the case where the observation points P; and P 4 * are at 
the same geometric distance a, the elements can be shown to 
depend only on a and the angular distance between P^ and Pj . 
This corresponds to the low low case with circular orbits. 
The positive definite matrix can be diagonalized with an 
orthogonal matrix with 




(42 


where 


A 


* diag(A,A ,A ,A , 

/ ' T 


' yyyyy > > > * ^ >0 (43 


The orthogonal matrix can be calculated by standard 
numerical techniques due to Jacobi. We can now define a set 
of orthogonal base functions. 


<l/.m 



p) 


(44 
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and write 


T (P) = 


where 


E a ^- (p) 


4- / 


3 * — 


1 E's t ( p. 


' s/ '\i 


(45 


(46 


The standard error of a^ is . Then the contribution 
of a* i pi to T can be explicitly determined, and only those 
base functions that are significant need be included. The 
choice of the set of base functions, however, is somewhat a 
matter of taste. Generally, the function becomes more 
oscillatory as i increases, i.e. for high wavenumbers. 
From this, we can use the size of to select the degree 
of smoothing. The economy of the inverse method hinges on 
being able to find an analytical expression for equation 
(41) , which will be described next. 

To be useful, the computation of should be fast 
enough so that an arbitarily large number of observations 
could be used. With n observations, or observation points, 
one has n(n+l)/2 values of to evaluate. Since is 
defined as a double integral over the sphere, evaluating 
this function could become very costly, using strictly 
numerical methods. Therefore we have begun by developing an 
algorithm to make this calculation simpler. The details are 
given in Lane and Gaposchkin (1986) . We want to compute 



■n^V t 


*7T P f * 


where 


-ir 


4*7T P? 


ft 5) v t 0 ©4A c40 ( 47 


> V ° & 


(48 


If 
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where 0;*^ is the spherical angle between ( 0 t * , X^ ) and 
y=( 0 , \ ) for arbitrary i f j. We may rewrite cos(0. ) as 


COO 0. 


- 603 0. . 0. r 

cj t) 


0*^.0. Ocej Q 


(49 


where a is the angle between the arcs 0^ and 0^ . Now as 
0 and X range over the whole sphere, a and 0. will 

also range over the whole sphere with 0$ a $ 2n * and 
O^0/^^ 7r » Hence, integrating the double kernel with 
respect to a and 0* will yield an equivalent integral, 
and we can write 



This can be rewritten as 



(51 


with 

b, * c, (. Cx «,Ct V'tJ.f)') 3 ’'*- 

c, - (r « v ) (. -r >- r x ) / i c -ft - 1 


Cl I <•>*-<** 

C>- - -t-TfH 

Cl - ^ 

C4. «• 

C*. * -t-fiH 

C» * ^ tli 

«*, 4 M = C.4. ‘•"J 4 

, <^W= Cj^C;,-C 4 


and 

V> t r ^ J 

= Z A* 
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(53 


So we can write the desired function as 

Jr 


n. 

u 


where 



bp) Y(ip) dip 


'Tf/jL 

Y(i P) = — / 4/1 +— OH sin 2 y 

H^vA-*** I f 


dy 


(54 


(55 


which can be written as an elliptic integral of the second 
kind as 


Y(0) 


_±_ 

14* 


E ( k ) 


(56 


where 


~Z = _ ^ ) 

1 - V?i( ip ) 


(57 


Rapidly convergent series for Elliptic Integrals are readily 
available. We have reduced the calculation of the double 
integral, to a single quadrature, which can be effectively 
done with standard proceedures such as a Romberg 
extrapolation. The accuracy of calculation of these matrix 
elements depends on the accuracy of the quadrature used, and 
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the user must balance speed with required accuracy. 

A further simplification and speed up is possible for 
the special case where z\, = r; . In this case the function 
depends only on the angular distance between the points. 
One could precompute these functions for a number of values 
of the separation and interpolate in this table. 


To investigate the effectiveness of this inverse method 
for the downward continuation problem, a sample calculation 
was performed. The details of the calculation are given in 
Lane, and Gaposchkin 1985) . The experiment consisted of 
selecting a grid of 24 mass points on a 1 deg x 1 deg grid 
at the earths surface. The mass points were selected to 
have zero average. They were arranged to have equal 
magnitude and alternate in sign. The center point of the 5 
x 5 grid was zero. The potential was computed, exactly, at 
121 grid points at an altitude of 212 km, on a 1 deg x 1 
deg, 11 x 11 grid. The potential was continued back down to 
an intermediate integration surface using this inverse 
method. The potential at this surface was recovered with an 
accuracy of 4%. 

In this experiment we had the advantage of knowing the 
true value of the potential by simply computing the 
potential from the mass points in the same way the observed 
data were computed. Hence, we were able to compute a 
measure of the performance of the inverse theory. For each 
test point (R,y ) we computed a kind of performance index, 
the standard deviation divided by the average value of the 
actual potential over the range of sample points 


s= 



C T (*) - 
E Ti \ 0 

4 . 


(58 


It is clear that the inverse method is computationally 
expensive, and so we give some information about the 
execution time of the routines. The computations were done 
of a Harris H800 computer, and all programming was done in 
FORTRAN. Table 1 lists the time involved for execution of 
principal steps for different values of N, the number of 



1 
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sample points. 

Table 1 . 

Computation time for Geopotential Inverse 

N CPU for CPU for Eigenvectors and 

Eigenvalues 


20 

18.471 sec 

0.417 

sec 

50 

113.708 sec 

6.309 

sec 

100 

449.639 sec 

48.136 

sec 

From table 

1 we can see 

that the 

computation 


eigenvalues and eigenvectors is slightly less than an n**3 
problem. This is consistent with the computation of the 
eigenvalues and eigenvectors being equivalent to computation 
of the inverse of . The computation of itself appears to be 
slightly less than an n**2 problem, but the computation of 
takes much longer. However, the computation time will 
change if one changes the sphere of integration and the 
sphere at which the sample points are taken, or if one 
changes the configuration of the sample points. For example 
the C£U time involved for computing where the 121 sample 
points are at 424 km height and the sphere of integration is 
212 km is only 13.523 sec: more than 30 times faster than 
the lower configuration! The computations in the remainder 
of the routine are linear, and so if one needs a 
satisfactory upper bound on the total CPU time involved for 
the 212 km height, 100 point configuration, then one can 
estimate it from (1/3) (n+n**2+n**3) 529 sec, where n=N/100. 

The numerical experiments conducted are described in 
detail in Lane and Gaposchkin (1986) . In partucular the 
spectral structure of the inversion. Therefore we briefly 
summarize the results in table 2. This test case concerns 
the 121 sample points on a sphere of r=l.lae, (height 
=637.8134 km) and and 25 test points on a sphere of R=1.05ae 
(height = 318.906) . 


Table 2 

Downward Continuation with 121 sample points 

Radius km S 

Sample Points 1.1 ae 637.813 0.019 

Test Points 1.05 ae 318.907 


Sample Points 1.033ae 212.604 


0.04 


II 
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Test Points 1.017ae 106.302 


Not only are these results very satisfactory, but 
of the program needed is much smaller than for 
computation in a global field. 


the size 
a similar 



If 
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GEOS-3/ATS-6 Tracking data analysis 

To test some of the algorithms for computing potential 
from SST data, the existing set of data from the GEOS-3 to 
ATS-6 tracking experiment was analysed. The daa was 
provided by J. Marsh (GSFC) and is described in Marsh et al 
1977. The data was provided in two files. For each of the 
observed 179 revolutions, Marsh computed an orbit based on 
all available tracking data, using the GEM10 gravity field, 
truncated at and complete through degree and order 12 as a 
reference. The residuals from SST observations were then 
computed. We interpret those residuals to be completely due 
to the anomalous potential above degree and order 12. Marsh 
then used a smoothing algorithm to obtain a smoothed range 
rate observation, and its rate or range acceleration. These 
three quantities, observed range rate, smoothed range rate, 
and smoothed range acceleration, were then given as a time 
tagged time series for each pass. In addition, for each 
pass, the state vectors for GEOS-3 and ATS-6 at a time near 
the beginning of the pass were given. The state vectors 
were presented as position and velocity in an inertial 
reference system. These state vectors allowed us to compute 
an approximate position for both satellites for each 
observation point by numerical integration. This provided 
the position and velocity used in evaluating the integrals 
for peturbations due to he anomalous potential. 

The GEOS-3 ATS-6 experiment was not ideal from our 
standpoint for a number of reasons. First, the data is not 
of 1 micron/sec limited by system parameters, and refraction 
effects. Second, the GEOS-3 satellite is at 800km altitude. 
Both of these factors limit the additional geopotential 
information to be gained from analysis of the data. In 
addition, the coverage was a fraction of one hemisphere. 
Also, the satellites were not surfaceforce compensated. 
Marsh took every possible step to account for 
nongravitational effects. Yet we believe there are still 
very important nongravitational effects in the reduced range 
rate residuals. Most important is, the High low 
configuration only measures a component of the along track 
velocity. When the lower satellite is just beneath the high 
satellite, it measures none. On the otherhand, the radial 
acceleration, which can be derived from numerical 
differentiation of the velocity data, directly measures the 
gravity acceleration. For the data taken at or near the sub 
satellite point, successful recovery of the gravity anomaly 
at satellite altitude has been demonstrated (Marsh et al, 
1977). However, for data away from the sub satellite point, 
a significant component of the along track velocity is 
detectable, and the question is, how to extract this from 
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the other velocity components. 

In any event, we proceed to set up the following 
analysis. Having solved the forward problem, in that, given 
an anomalous potental T, we can compute the velocity in all 
components, these can be projected on the line of sight from 
ATS-6 to GEOS-3 . To obtain, a correction to the anomalous 
potential, we assume that the observed residual velocity on 
any iteration, is all along track, and is therefore a 
measure of the potential at that point in space. 
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Sununary and Conclusions 

A relatively efficient computation has been developed 
for the kernal function used in the Geophysical Inverse 
Problem. 

Geophysical Inverse Theory seems a viable method to 
solve the downward continuation problem, and map the 
geopotential from satellite altitudes to the Earth's 
surface. 

Geophysical Inverse Theory seems able to recover the 
geopotential at the Earth's surface with an accuracy of a 
few percent, in the absence of noise. 

The Inverse problem seems appropriately broken up into 
local networks, which can be solved separately, leading to 
significantly simplification in computing requirements. 


The analytical development for computing perturbations 
in velocity have been developed and are shown to be have 
sufficient accuracy for the analysis of SST data. 


Future Studies 

For the Geophysical Inverse Method 

The trade-off of size of the sample point space and the 
solution space. 

The influence of errors in the sample point space, both 
position and measurement errors. 

The benefit of sample points at different altitudes, in 
improving the solution, and eliminating possible 
nonuniqueness. 

Formal understanding of the annihilator, and possible 
tests for nonuniqueness. 


Most important, it is essential to complete this 
analysis by merging the Satellite perturbation analysis with 
the downward continuation calculation, and to iterate on 
this process until convergence. Each part is demonstrated 
to be effective. 
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